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Abstract 



We consider the numerical solution of the equation — Au— /(«) = g, for the 
' unknown u satisfying Dirichlet conditions in a bounded domain f2. The non- 

, linearity / has bounded, continuous derivative. The algorithm uses the finite 

element method combined with a global Lyapunov-Schmidt decomposition. 
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> ! 1 Introduction 

O i 

| We consider the partial differential equation 
^ 1 F(u) = -Au - f(u) = g, u\ dn = 0, 

m ' 

on domains il e R", taken to be open, bounded, connected subsets of M. n with 
piecewise smooth boundary <9f2, assumed to be at least Lipschitz at all points. 
There is a vast literature concerning the number of solutions for general and positive 
solutions for different kinds of nonlinearity / and right-hand side g (to cite a few, 

m, h, [a, m, 0, i, h, i, i). 

| Here we assume that the nonlinearity / : R — » K has a bounded, continuous, 

derivative, a < f'(y) < We show how a global Lyapunov-Schmidt decomposi- 
tion introduced by Berger and Podolak [TU] in their proof of the Ambrosetti-Prodi 
theorem ([3], [1]) gives rise to a satisfactory solution algorithm using the finite ele- 
ment method. The decomposition was rediscovered by Smiley , who realized its 
potential for numerics: our results advance along these lines. 

Write Ad for the Dirichlet Laplacian in fi. The algorithm is especially conve- 
nient when the number d of eigenvalues of — Ab in the range of /' is small: the 
infinite dimensional equation reduces to the inversion of a map from M. d to itself. 

The subject of semilinear elliptic equations is sufficiently mature that algorithms 
should stand side by side with theory. The situation may be compared to the study 
of functions of one variable in a basic calculus course. Some functions, like parabo- 
las, may be handled without substantial computational effort, but understanding 
increases with graphs, which are obtained by following a standard procedure. 

We do not handle the difficulties and opportunities related to the finite dimen- 
sional inversion: a generic solver (as in [12] and, for d = 2, [13]) should be replaced 
by an algorithm which makes use of features inherited by the original map F. Here 
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we only deal with examples for which d = 1 and 2, and there is some craftsmanship 
in handling the 2-dimensional example. It is in this step of the PDE solver that 
delicate issues like nonresonance and lack of properness come up. 

The results are part of the PhD thesis of the first author [TJj . Complete proofs 
are presented elsewhere. The authors are grateful to CAPES, CNPq and Faperj for 
support. 

2 The basic estimate 

We consider the semilinear elliptic equation presented in the introduction for a 
nonlinearity / : R — > K with bounded, continuous derivative. 

With these hypotheses, it is not hard to see that F(u) = —Au — f(u) is a C 1 
map between the Sobolev spaces Hq(£1) and L 2 (fl) — H°(£l) and between H^fl) 
and i/ _1 (f2) ~ Hq(H). We concentrate on the second scenario, which is natural for 
the weak formulation of the problem. Still, the geometric statements below hold in 
both cases. To fix notation, set F : X — >• Y, where X — H^(ft) and Y = H~ 1 (fl). 

The basic estimate is given in Proposition [TJ Its proof is a simple extension of 
the argument in |10j . 

Define /'(R) = [a, b] (a allowed to be -co) and a larger interval [5,6] D [a, b]. 
Label the eigenvalues of —Ad in non-decreasing order. The index set J associated 
to [a,b] is the collection of indices of eigenvalues of —Ad in that interval. The set 
J is associated to the nonlinearity / if [5, b] = /'(R). An index set defined this way 
is complete: it contains all indices labeling an eigenvalue in the interval. 

Denote the vertical subspaces by V x C X and V Y C Y the spans of the normal- 
ized cigenfunctions 4>j,j € J in X and Y respectively, with orthogonal complements 
W x and W Y . Let P and Q be the orthogonal projections on V and W. Clearly, the 
dimension of the vertical subspaces equals \J\, the cardinality of J. Let v + W x C X 
be the horizontal affine subspace of vectors v + w, w € W x and consider a projected 
restriction F v : v + W x — > W Y , the restriction of P Y F to v + W x . 

Proposition 1. Let J be the index set associated to the nonlinearity f (or to any 
interval [5,6] containing f'(M)). Then the derivatives DF V : v + W x — > W Y are 
uniformly bounded from below. More precisely, there exists C > such that 

Vug V x Vwev + W x VheW x , \DF v {w)h\ Y >C\h\ x . (1) 

All such maps are invertible. 

A direct application of Hadamard globalization theorem ([15]) implies that the 
projected restrictions are diffeomorphisms, for each v G V x . 

3 The underlying picture 

The geometric implications are very natural. The image under F of each horizontal 
affine subspace v + W x is a sheet, i.e., a surface which projects under P Y diffcomor- 
phically to the horizontal subspace W Y . In particular, every vertical affine subspace 
w + V Y intercepts each sheet exactly at a single point. It is not hard to see that 
the intersection is transversal: tangent spaces of sheet and vertical affine subspace 
form a direct sum decomposition of Y. 

A fiber is the inverse image of a vertical affine subspace. In a similar fashion, 
fibers are surfaces of dimension \J\ which meet every horizontal affine subspace 
v + W x at a single point — again, the intersection is transversal. Thus, a vertical 
subspace parameterizes diffeomorphically each fiber, or, said differently, each fiber 
has a single point of a given height. 
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Recall a key idea in [TU] and [12] ■ It is clear that X and Y are respectively 
foliated by fibers and vertical affine subspaces. By definition, all the solutions of 
F(u) — g must lie in the fiber a g — F^ 1 (g + V Y ). So, in principle, one might solve 
the equation by first identifying a g ~ M' 17 ' and then facing the finite dimensional 
inversion of F : a g — > g + V Y ■ 

Horizontal affine subspaces are taken diffeomorphically to sheets, but fibers are 
not taken diffeomorphically to vertical affine subspaces. In a sense, the nonlinearity 
of the problem was reduced to a finite dimensional issue. 



Figure 1: Horizontal affine subspace, fiber; sheet, vertical affine subspace 



4 Finding the fiber 

Recall that each horizontal affine subspace v + W x contains exactly one element of 
each fiber. So, to identify a g , choose v + W x and search in it for an element of a g . 
Said differently, one may think of F v : v + W x —> W Y as being a diffeomorphism 
between fibers (represented by points in v + W x ) and vertical affine subspaces 
(represented by points in W Y ). The situation is ideal for an application of Newton's 
method: local improvements are performed by linearization of the diffeomorphism. 

There is one difficulty, however, related to implementation issues. The functional 
spaces X and Y give rise to finite dimensional vector spaces generated by finite 
elements. We provide some detail: an excellent reference is [IB] . First of all, 
triangulate the domain Q, i.e., split it into disjoint simplices in R™. In the examples 
in Section[71 fl is the uniformly triangulated rectangle [0, 1] x [0, 2]. A nodal function 
is a continuous function that is affine linear on each simplex and has value one at a 
given vertex and zero at the remaining vertices. These functions form a nodal basis, 
which spans a finite dimensional subspace of i/p(Sl). Figure [2] shows an example of 
a triangulation of f2 and one nodal function. 

Inner products of nodal functions, both in H 1 and L 2 are often zero, a fact which 
simplifies the numerics associated to the weak formulation of the equation F(u) = g. 
The vertical subspaces V x and V Y are spanned by eigenfunctions <pj ,j € J and are 
well approximated by a few linear combinations on the nodal basis. On the other 
hand, obtaining a similar basis for the approximation of the orthogonal subspaces 
W x and W Y requires much more numerical effort and should be avoided. 

To circumvent this problem, extend the Jacobian of F v : v + W x — > W Y at a 
point u to an invertible operator L u : X — > Y which is easy to handle and apply 
Newton's method to L u instead. Setting 



it is clear that L u has the required properties: it takes W x to W Y and V x to V Y and 
the restriction to v + W x equals DF V , which is invertible. Moreover, the restriction 




L u z = -Az - P Y f'{u)P x z, 
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Figure 2: Uniform triangulations of [0, 1] x [0, 2] and a nodal function 

to V x coincides with —A. This map is no longer a differential operator, due to the 
integrals needed to compute the projections P. But those new terms are innocu- 
ous in the finite element formulation — the sparsity of the underlying matrices is 
preserved, together with the possibility of standard preconditioning routines. 




Figure 3: Finding the right fiber 

We search for a point of a horizontal affine subspace v + W x which belongs to 
a 9 ,g £ y.The algorithm is straightforward: see Figured] Choose a starting point 
Mo and consider its image Fo. All would be well if the projections of Fq and g 
on the horizontal subspace W Y were equal or at least very near. When this does 
not happen, proceed by a continuation method to join both projections. Notice 
that the algorithm searches for the fiber (i.e., for a point in the fiber) by moving 
horizontally in the domain. A direct Newton iteration does not work necessarily: 
think of finding the (trivial) root of arctan(x) = starting sufficiently far from the 
origin. 

5 Moving along the fiber 

The necessary ingredients for a simple predictor-corrector method to move along a 
fiber are now available. Say u g a g and we want to find another point in a g . Recall 
that fibers are parameterized by height v e V x . Take u + v, which is probably not 
in a g , as a starting point for the algorithm in Section [4] to obtain the point of a g 
in the same horizontal affine subspace of u + v (see Figure 0] for two such steps). 

We don't know much about the behavior of F restricted to a fiber: the hypothe- 
ses on the nonlinearity / are not sufficient to imply properness of F, for example. 
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Figure 4: Mapping a 1-D fiber 



In particular, it is not clear that the restrictions of F to fibers are also proper. 

6 Stability issues 

Proposition Q] in Section [2] ensures geometric stability, in the sense that the global 
Lyapunov-Schmidt decompositions preserve their properties under perturbations. 
This is convenient when replacing the vertical subspaces spanned by eigenfunctions 
by their finite clement counterparts. 

As for the algorithm itself, the identification of the fiber is robust, being a 
standard continuation method associated to a diffcomorphism between horizontal 
affine spaces. The numerical analysis along a fiber is a different matter, and the 
fundamental issue was addressed by Smiley and Chun 12 : they showed that the 
finite element approximations to the restriction of the function F to (compact sets 
of) the fiber can be made arbitrarily close to the original map in the appropriate 
Sobolev norm. Here one must proceed with caution: a small metric perturbation 
may alter the number of solutions, as when switching from x i— > x 2 ,x G R to 
x H» x 2 — e, which is a perturbation of order e for arbitrary C k norms. Still, solutions 
of F which are regular points are stable: they correspond to nearby solutions of 
sufficiently good finite element approximations F . 

A different approach might be to interpret the algorithm as a provider of good 
starting points for Newton's iteration or at least a continuation method. As stated 
in |17j . computer assisted arguments require good approximations for the eventual 
validation of solutions. 

7 Some examples 

For the examples that follow, F(u) = — Au — f(u) = g, with Dirichlet conditions 
on 17 = [0,1] x [0,2]. Here, —Ad has simple eigenvalues and Ai = jir 2 « 12.34, 
A 2 = 2ir 2 ss 19.74, A 3 = ^-tt 2 » 41.95. Denote by 4>* and cf)\ the eigenfunctions of 
—Au normalized in X and Y . 

The first example is a nonlinearity / satisfying the hypotheses of the Ambrosetti- 
Prodi theorem with /(0) = and derivative fix) = aarctan(a;) + f3 with 



The graph of /' is shown on the left of Figure [5j Here, the index set associated 
to / is J = {!}: V x and V Y are spanned by <j)f ,(f>Y > 0. For right-hand side set 
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Figure 5: The derivative of / and the image of a g 

g{x) — — 100x(:r — 1)2/(2/ ~ 2), shown on the left of Figure [51 which has a large 
negative component along the ground state. 

We search for an element of the fiber a g in the horizontal subspace W x , starting 
from uq = 0, in the notation of Section @] The result is the function on the right of 
Figure [BJ Now move along a g , as in Section [5j The graph on the right of Figure [5] 




Figure 6: A right-hand side and a function on its fiber 

plots (u,(/)*) x (the height of u 6 a g ) versus (F(u), <f)\ ) Y (the height of F{u)). 
The horizontal line indicates the height of g: the solutions of the original PDE 
correspond to the intersections between the curve and this line. The two solutions 
found in this case are presented in Figure [7l 




Figure 7: Ambrosetti-Prodi Solutions 
For the next example, J = {1} but / is a nonconvex function whose derivative is 



G 



depicted on the left of Figure[5] We consider the fiber through Uo(x) = —50 <j>x{x) + 
10 4>2{x), i.e., a F f Uo y According to Figure[5] (right), moving up the fiber yields three 
distinct solutions. 




Figure 8: Non-convex / 



As a concluding example, we take a nonlinearity / for which J = {1,2}: here 
the vertical spaces are spanned by the first two eigenfunctions. The function / is 
of the same form as the first example and its derivative is shown in Figure [9] We 
study the fiber through the point uq = 0, which is cto, since F(0) = 0. Recall from 
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Figure 9: The range of /' contains Ai and A2 

Section |4] that there is exactly one point of ao for each height, i.e., given a point 
u € V x , there is a unique point £(u) £ ao in the same horizontal affine subspace as 
u. For a circle C centered at the origin in V x , ({C) <E ao. The image F(((C)) is 
shown in the right side of Figure llOl here, we must project F(£(u)) along directions 
<pY and 02 '• Clearly, there is a double point Z in F(((C)) and it is not hard to 
identify in C its two pre-images, U and D marked in the left of Figure [TU] 
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Figure 10: Two preimages on the circle 

We now obtain two additional preimages of Z in a rather naive fashion. The 
images under F o ( of the four half-axes of V x are drawn on the right of Figure [HI 
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Figure 11: Two preimages along the horizontal axis 




Figure 12: The four solutions 
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